On part du Lagrangien classique $$L=\frac{m}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2)-mgz,$$
Puis on introduit les coordonnée sphériques $$ \begin{cases} x = r\sin(\theta)\cos(\phi)\\ y = r\sin(\theta)\sin(\phi)\\ z = r\cos(\theta) \end{cases} $$ où $r$ est constant mais $\theta$ et $\phi$ dépendent du temps.
On arrive alors au Lagrangien $$ L=\frac{mr^2}{2}(\dot{\theta}^2+\sin^2(\theta)\dot{\phi}^2)-mgr\cos(\theta). $$
Les équations d'Euler-Lagrange sont maintenant $$ \begin{cases} \ddot{\theta} = \sin(\theta)(\cos(\theta)\dot{\phi}^2+\frac{g}{r})\\ \sin^2(\theta)\left(\ddot{\phi}+\frac{2\dot{\theta}\dot{\phi}}{\tan(\theta)}\right) = 0 \end{cases} $$
On introduit les coordonnées $$ \begin{cases} \alpha = mr^2\dot{\theta}\\ \beta = mr^2\sin^2(\theta)\dot{\phi} \end{cases} $$
Le Hamiltonien est alors $$ H:=\dot{\theta}\alpha+\dot{\phi}\beta-L=\frac{\alpha^2}{2mr^2}+\frac{\beta^2}{2mr^2\sin^2(\theta)}+mgr\cos(\theta) $$
Finalement on peut écrire les équations de Hamilton $$ \begin{cases} \dot{\theta} = \frac{\partial H}{\partial \alpha } = \frac{\alpha}{mr^2}\\ \dot{\phi} = \frac{\partial H}{\partial \beta } = \frac{\beta}{mr^2\sin^2(\theta)}\\ \dot{\alpha} = -\frac{\partial H}{\partial \theta} = mgr\sin(\theta)+\frac{\beta^2 \cos(\theta)}{mr^2\sin^3(\theta)}\\ \dot{\beta} = -\frac{\partial H}{\partial \phi } = 0 \end{cases} $$
On a une singularité dans le Lagrangien dûe au $\sin^2(\theta)$ au dénominateur.
Cependant comme $H$ est une constante du mouvement, en prenant compte de la positivité des deux premiers termes et du caractère borné du troisième on n'est pas censé avoir explosion du point de vue théorique et numériquement tant que le Hamiltonien est à peu près conservé non plus.
import numpy as np
from scipy.integrate import odeint
m, g, r = 1, 10, 1
def flux(X, t):
theta, phi, alpha, beta = X
dtheta = alpha/(m*r**2)
dphi = beta/(m*r**2*(np.sin(theta))**2)
dalpha = m*g*r*np.sin(theta)+beta**2*np.cos(theta)/(m*r**2*(np.sin(theta))**3)
dbeta = 0
return np.array([dtheta, dphi, dalpha, dbeta])
ini = np.array([np.pi/2, 0, 0.5, 0.5])
temps = np.arange(0, 50, 0.05)
resultat = odeint(func=flux, y0=ini, t=temps)
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import HTML
%matplotlib inline
plt.rc("figure", figsize=(12,9))
plt.rcParams["animation.writer"] = "avconv"
fig = plt.figure()
ax = Axes3D(fig)
theta, phi =resultat[:,0], resultat[:,1]
x = r*np.sin(theta)*np.cos(phi)
y = r*np.sin(theta)*np.sin(phi)
z = r*np.cos(theta)
ligne, = ax.plot(x, y, z, ls="--", color="blue")
fig
fig = plt.figure()
ax = Axes3D(fig)
traj, = ax.plot3D([], [], [], ls='--', color="blue")
tige, = ax.plot3D([], [], [], color="red", marker="o")
ax.set_xlim([np.min(x), np.max(x)])
ax.set_ylim([np.min(y), np.max(y)])
ax.set_zlim([np.min(z), np.max(z)])
ax.axes.set_aspect("equal")
ax.view_init(20, 60)
def maj(k):
traj.set_data(x[:k], y[:k])
traj.set_3d_properties(z[:k])
tige.set_data([0, x[k]], [0, y[k]])
tige.set_3d_properties([0, z[k]])
ax.view_init(20, 60-0.1*k)
return traj, tige
ani = FuncAnimation(fig=fig, func=maj, frames=range(len(x)), interval=50, blit=True)
HTML(ani.to_html5_video())